import numpy as np
import pandas as pd
import scipy.signal as signal
import matplotlib.pyplot as plt
import pywt

# 生成电流有效值波形
def current_to_I_rms(I):
    # 窗口大小和步长
    window_size = 0.02  # 窗口大小，单位：秒
    step_size = 0.005  # 步长，单位：秒
    # 计算窗口数量
    num_windows = int((signal_duration - window_size) / step_size) + 1
    # 初始化有效值数组
    I_rms_over_time = np.zeros(num_windows)
    # 滑动窗口计算有效值
    for i in range(num_windows):
        start_index = int(i * step_size * sampling_frequency)
        end_index = int((i * step_size + window_size) * sampling_frequency)
        window = I[start_index:end_index]
        I_rms_over_time[i] = np.sqrt(np.mean(window ** 2))
    # 计算每个窗口的中心时间点
    window_times = np.linspace(0, signal_duration - window_size, num=num_windows)
    return I_rms_over_time, window_times
#对有效值波形进行过滤
def filter_of_I_rms(I_rms_over_time):
    fs = 200  # 采样频率
    cutoff = 2  # 截止频率
    order = 4  # 滤波器阶数
    b, a = signal.butter(order, cutoff, fs=fs, btype='low')

    filtered_I_rms = signal.filtfilt(b, a, I_rms_over_time)
    return filtered_I_rms
def wavelet_transform(wavelet_data, wavelet, level):
    coeffs = pywt.wavedec(wavelet_data, wavelet, level=level)
    reconstructed_signal = pywt.waverec(coeffs, wavelet)
    return coeffs, reconstructed_signal

def select_wavelet_threshold(coeffs):
    leng = len(coeffs[1])
    print(leng)
    I_rms_diff = np.diff(coeffs[0])
    # 阈值设置，在此选择功率变换最小的灯泡为阈值
    # threshold = 1.4 * np.std(coeffs[1])
    lower_limit = np.percentile(I_rms_diff, 98)
    upper_limit = np.percentile(I_rms_diff, 2)
    # threshold = max(abs(lower_limit), abs(upper_limit))
    threshold = 2
    # 找出超过阈值的系数点的索引
    t_shijian = np.arange(0, signal_duration, signal_duration / leng)
    exceeding_indices = np.where(np.abs(I_rms_diff) > threshold)[0]

    # 确定时间点间隔阈值，例如可以是时间点数组的均值间距的某个百分比
    interval_threshold = 4

    # 初始化代表性时间点列表和临时时间点列表
    representative_points = []
    temp_points = []

    # 遍历时间点，分组连续的时间点
    for i in range(len(exceeding_indices)):
        if len(temp_points) == 0 or (exceeding_indices[i] - temp_points[-1]) <= interval_threshold:
            # 如果是第一个时间点或当前时间点与最后一个时间点间隔小于阈值，则加入临时列表
            temp_points.append(exceeding_indices[i])
        else:
            # 否则，将临时列表的时间点加入代表性时间点列表，并重置临时列表
            representative_points.append(temp_points[0])  # 添加中间时间点或第一个时间点
            temp_points = [exceeding_indices[i]]

    # 添加最后一组时间点
    if temp_points:
        representative_points.append(temp_points[0])

    # selected_time_points = t_shijian[representative_points]
    selected_time_points = [i / leng * signal_duration for i in representative_points]
    data_points = [int(i * 3200) for i in selected_time_points]
    print(exceeding_indices)
    print("Selected time points with significant detail coefficients:", selected_time_points)
    print(representative_points)
    print(data_points)
    return selected_time_points, data_points

# 参数设置
sampling_frequency = 3200  # 采样频率，单位：Hz
signal_duration = 115         # 信号时长，单位：秒

# 生成时间向量
name_of_txt = 'current20'
# index0 = 2
t = np.arange(0, signal_duration, 1/sampling_frequency)
df = pd.read_csv(f'{name_of_txt}.csv', encoding='utf-8')

# 假设我们要读取名为'ColumnName'的列
column_name = 'after'
column_data = df[column_name].tolist()  # 将列数据转换为列表
# 生成电流信号，这里我们使用一个简单的正弦波作为示例
endindex=signal_duration*3200
I = np.array(column_data[0:endindex])  # 50Hz的基波

I_rms_over_time, window_times = current_to_I_rms(I)
filtered_I_rms = filter_of_I_rms(I_rms_over_time)

wavelet = 'db4'
level = 3
coeffs, reconstructed_signal = wavelet_transform(filtered_I_rms, wavelet, level)
# 从系数重建信号

# 计算D4的一阶差分
cD4_diff = np.diff(coeffs[1])
I_rms_diff = np.diff(coeffs[0])


# df = pd.DataFrame(I_rms_diff, columns=['Detail Coefficients D4'])
# df.to_csv(f'detail_coefficients_{name_of_txt}.csv', index=False)

# 事件检测函数，通过小波变换后的近似系数的一阶差分，设置某个阈值进行挑选
selected_time_points, data_points = select_wavelet_threshold(coeffs)

start_index1 = data_points[index0]-3200
start_index2 = data_points[index0]+3200
temp_signal_duration = 0.1
length = int(temp_signal_duration*sampling_frequency)

current1 = np.array(column_data[start_index1:start_index1 + length])
current2 = np.array(column_data[start_index2:start_index2 + length])
current3 = current2 - current1
t_temp = np.arange(0, temp_signal_duration, 1/sampling_frequency)

# 计算傅里叶变换
fft_current1 = np.fft.fft(current1)
frequencies1 = np.fft.fftfreq(len(current1), d=1/sampling_frequency)

fft_current2 = np.fft.fft(current2)
frequencies2 = np.fft.fftfreq(len(current2), d=1/sampling_frequency)

fft_current3 = np.fft.fft(current3)
frequencies3 = np.fft.fftfreq(len(current3), d=1/sampling_frequency)

#运行相位补偿法，计算先添加的负载波形
phase_diff = np.angle(fft_current2) - np.angle(fft_current1)
fft_before_avg_corrected = np.abs(fft_current1) * np.exp(1j * (np.angle(fft_current1) + phase_diff))
fft_new_load = fft_current2 - fft_before_avg_corrected
new_load_waveform = np.fft.ifft(fft_new_load)
new_load_waveform = np.real(new_load_waveform)
# 绘制原始电流波形和随时间变化的有效值波形
def plot():
    plt.figure(figsize=(12, 12))

    plt.subplot(4, 2, 1)
    plt.plot(t, I, label='Current Waveform')
    plt.title('Current Waveform')
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude (A)')
    plt.legend()

    plt.subplot(4, 2, 2)
    plt.plot(window_times, I_rms_over_time, label='RMS over Time', color='red')
    plt.title('RMS Value Over Time')
    plt.xlabel('Time (s)')
    plt.ylabel('RMS Amplitude (A)')
    plt.legend()

    # plt.subplot(4, 2, 3)
    # plt.plot( filtered_I_rms, label='Filtered Signal (Moving Average)')
    # plt.legend()

    # 绘制第一级近似系数
    plt.subplot(4, 2, 3)
    plt.plot(coeffs[0])
    plt.title(f'Level {level} Approximation Coefficients')
    plt.gray()

    # 绘制第一级细节系数
    plt.subplot(4, 2, 4)
    plt.plot(coeffs[1])
    plt.title(f'Level {level} Detail Coefficients')
    plt.gray()
    #
    plt.subplot(4, 2, 5)
    plt.plot(t_temp, current3, label='Current Waveform')
    plt.title('Current Waveform')
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude (A)')
    plt.legend()

    plt.subplot(4, 2, 6)
    plt.plot(t_temp, new_load_waveform, label='Current Waveform')
    plt.title('Current Waveform')
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude (A)')
    plt.legend()

    # 绘制第一级近似系数
    plt.subplot(4, 2, 7)
    plt.plot(I_rms_diff)
    plt.title(f'Level {level - 3} Detail Coefficients')
    plt.gray()

    # 绘制第一级细节系数
    plt.subplot(4, 2, 8)
    plt.plot(cD4_diff)
    plt.title('Level 3 Detail Coefficients')
    plt.gray()

    plt.tight_layout()
    plt.show()

    # # 计算需要的子图数量
    # num_subplots = 1 + len(coeffs) * 2
    # plt.figure(figsize=(15, 10))
    #
    # # 绘制原始信号
    # plt.subplot(9, 1, 1)
    # plt.plot(t, I)
    # plt.title('Original Signal')
    #
    # print(num_subplots)
    # # 绘制近似系数和细节系数
    # for i, (a, d) in enumerate(zip(coeffs[::2], coeffs[1::2]), 1):
    #     plt.subplot(9, 2, 2 * i - 1)
    #     plt.plot(a)
    #     plt.title(f'Approximation Coefficients at Level {i}')
    #
    #     plt.subplot(9, 2, 2 * i)
    #     plt.plot(d)
    #     plt.title(f'Detail Coefficients at Level {i}')
    #
    # # 显示图表
    # plt.tight_layout()
    # plt.show()

plot()